home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / oper_sys / weyl / weyl_lsp.lha / bigfloat.lisp < prev    next >
Encoding:
Text File  |  1991-10-04  |  40.8 KB  |  1,229 lines
  1. ;;; -*- Mode:Lisp; Package:Weyli; Base:10; Lowercase:T; Syntax:Common-Lisp -*-
  2.  
  3. ;; $Id: bigfloat.lisp,v 2.13 1991/10/04 22:42:51 rz Exp $
  4.  
  5. ;; Arbitrary Precision Real Arithmetic System
  6. ;;  by Tateaki  Sasaki
  7. ;;  The UNIVERSITY OF UTAH,  MARCH 1979
  8.  
  9. ;; For design philosophy and characteristics of this system, see T. Sasaki, 
  10. ;; ``An Arbitrary Precision Real Arithmetic Package In Reduce,''
  11. ;;  Proceedings of Eurosam '79, Marseille (FRANCE), June 1979.
  12.  
  13. ;; For Implementation notes and using this system, see T. Sasaki, 
  14. ;; ``Manual For Arbitrary Precision Real Arithmetic System in Reduce,' 
  15. ;; Operating Report of Utah Symbolic Computation Group
  16.  
  17. ;;  In order to speed up this system, you have only to rewrite four
  18. ;; routines (DECPREC!, INCPREC!, PRECI!, AND ROUND!LAST)
  19. ;; machine-dependently.
  20.  
  21. ;;  This function constructs an internal representation of a number
  22. ;; 'n' composed of the mantissa MT and the exponent EP with the base
  23. ;; 10.  The magnitude of the number thus constructed is hence
  24. ;; MT*10^EP.
  25.  
  26. ;; **** CAUTION!  MT and EP are integers.  So, EP denotes the order of
  27. ;;                the last figure in 'N', WHERE order(N)=K if 10**K <=
  28. ;;                ABS(N) < 10**(K+1), with the exception order(0)=0. 
  29. ;;
  30.  
  31. ;; The number 'n' is said to be of precision 'k' if its mantissa is a
  32. ;; k-figure number.  MT and EP are any integers (positive or
  33. ;; negative).  So, you can handle any big or small numbers.  in this
  34. ;; sense, 'BF' denotes a big-floating-point number.  Hereafter, an
  35. ;; internal representation of a number constructed by MAKE-BIGFLOAT is
  36. ;; referred to as a big-float representation.
  37.  
  38. (in-package "WEYLI")
  39.  
  40. (defmethod print-object ((d real-numbers) stream)
  41.   #+Genera
  42.   (format stream "~'bR~(~D)" (fp-precision d))
  43.   #-Genera
  44.   (format stream "R(~D)" (fp-precision d)))
  45.  
  46.  
  47. (defsubst make-bigfloat (domain mantissa exponent)
  48.   (make-instance 'bigfloat :domain domain
  49.          :mantissa mantissa :exponent exponent))
  50.  
  51. ;; This function returns t if x is a big-float representation, else it
  52. ;; returns NIL.
  53. (defun bigfloatp (x)
  54.   (eql (class-name (class-of x)) 'bigfloat))
  55.  
  56. (defmethod number? ((element bigfloat))
  57.   t)
  58.  
  59. ;; This function counts the precision of a bigfloat 'n'.
  60. ;;
  61. ;; FIXTHIS:  The 1+ below is total kludge.  This whole package needs
  62. ;; to be rewritten sometime soon.  
  63. (defun preci! (nmbr)
  64.   (let* ((mt (1+ (abs (bigfloat-mantissa nmbr))))
  65.      (len (integer-length mt)))
  66.     (if (< len 10)
  67.     (lisp:ceiling (lisp:log mt 10))
  68.     (lisp:ceiling (lisp:+ (lisp:/ (lisp:+ len -10) #.(lisp:log 10 2))
  69.                   (lisp:log (lisp:ash mt (lisp:- 10 len)) 10))))))
  70.  
  71. ;; This function counts the order of a bigfloat 'n'.
  72. ;; **** ORDER(N)=K IF 10**K <= ABS(N) < 10**(K+1) 
  73. ;; ****     WHEN N IS NOT 0, AND ORDER(0)=0.      
  74. (defun order! (nmbr)
  75.   (if (lisp:zerop (bigfloat-mantissa nmbr)) 0
  76.       (lisp:+ (preci! nmbr) (bigfloat-exponent nmbr) -1)))
  77.  
  78. (defun convert-number->characters (number)
  79.   (if (lisp:zerop number) (list #\0)
  80.       (let ((chars ()))
  81.     (loop with n = number
  82.           and digit
  83.           while (not (lisp:zerop n))
  84.           do (multiple-value-setq (n digit) (lisp:truncate n 10))
  85.          (push (digit-char digit) chars))
  86.     chars)))
  87.  
  88. (defmethod print-object ((number bigfloat) stream)
  89.   (setq number (round!mt number (lisp:- (fp-precision (domain-of number)) 2)))
  90.   (with-slots (mantissa exponent) number 
  91.     (let ((u (convert-number->characters (abs mantissa)))
  92.       k)
  93.       (flet ((bfprin1 (u)
  94.            (when (lisp:minusp mantissa)
  95.          (push #\- u))
  96.            ;; Suppress trailing zeroes
  97.            (loop for v on (reverse u)
  98.              while (and (char= (first v) #\0)
  99.                 (if (rest v) (char/= (second v) #\.)
  100.                     t))
  101.              finally (setq u (nreverse v)))
  102.            ;; Now print the number
  103.            (loop for char in u
  104.              do (write-char char stream))))
  105.     (or (rest u) (push #\0 (rest u)))
  106.     (push #\. (rest u))
  107.     (bfprin1 u)
  108.     (princ "B" stream)
  109.     (setq u (convert-number->characters (abs (setq k (order! number)))))
  110.     (setq u (cons (if (lisp:< k 0) #\- #\+)
  111.               u))
  112.     (loop for char in u
  113.           do (write-char char stream))
  114.     number))))
  115.  
  116.  
  117. ;;; Routines For Converting A Big-Float Number
  118.  
  119. ;; This function converts a number N to an equivalent number the
  120. ;; precision of which is decreased by K.
  121. (defun decprec! (number k)
  122.   (with-slots (exponent mantissa) number
  123.     (make-bigfloat *domain*
  124.            (lisp:truncate mantissa (lisp:expt 10 k))
  125.            (lisp:+ exponent k))))
  126.  
  127. ;;  This function converts a number N to an equivalent 
  128. ;;       number the precision of which is increased by K.
  129. (defun incprec! (number k)
  130.   (with-slots (exponent mantissa) number
  131.     (make-bigfloat *domain*
  132.            (lisp:* mantissa (lisp:expt 10 k))
  133.            (lisp:- exponent k))))
  134.  
  135. ;; This function converts a number to an equivalent number of precision
  136. ;; k by rounding or adding '0's.
  137. (defun conv!mt (number k)
  138.   (unless (bigfloatp number)
  139.     (error "Expected ~S to be a bigfloat" number))
  140.   (unless (and (integerp k) (lisp:> k 0))
  141.     (error "Expected ~S to be a positive integer" k))
  142.   (cond ((lisp:zerop (setq k (lisp:- (preci! number) k)))
  143.      number)
  144.     ((lisp:< k 0) (incprec! number (lisp:- k)))
  145.     (t (round!last (decprec! number (1- k))))))
  146.  
  147. ;; This function converts a number 'n' to an equivalent number having
  148. ;; the exponent k by rounding 'n' or adding '0's to 'n'.
  149. (defun conv!ep (nmbr k)
  150.   (unless (bigfloatp nmbr) 
  151.     (error "Invalid argument to conv!ep: ~S" nmbr))
  152.   (unless (integerp k)
  153.     (error "Invalid second argument to conv!ep: ~S" k))
  154.   (cond ((lisp:zerop (setq k (lisp:- k (bigfloat-exponent nmbr))))
  155.      nmbr)
  156.     ((lisp:< k 0) (incprec! nmbr (lisp:- k)))
  157.     (t (round!last (decprec! nmbr (lisp:- k 1))))))
  158.  
  159. ;; This function returns a given number N unchanged if its precision
  160. ;; is not greater than K, else it cuts off its mantissa at the (K+1)th
  161. ;; place and returns an equivalent number of precision K.
  162. ;;  **** CAUTION!  NO ROUNDING IS MADE.                  
  163. (defun cut!mt (nmbr k)
  164.   (declare (fixnum k))
  165.   (unless (bigfloatp nmbr)
  166.     (error "Invalid argument to cut!mt: ~S" nmbr))
  167.   (unless (and (integerp k) (lisp:> k 0))
  168.     (error "Invalid precision to cut!mt: ~D" k))  
  169.   (if (lisp:plusp (setq k (lisp:- (preci! nmbr) k)))
  170.       (decprec! nmbr k)
  171.       nmbr)) 
  172.  
  173. ;; This function returns a given number N unchanged if its exponent is
  174. ;; not less than K, else it cuts off its mantissa and returns an
  175. ;; equivalent number of exponent K.
  176. ;;  **** CAUTION!  NO ROUNDING IS MADE.                  
  177. (defun cut!ep (nmbr k)
  178.   (unless (bigfloatp nmbr)
  179.     (error "Invalid argument to cut!ep: ~S" nmbr))
  180.   (unless (integerp k)
  181.     (error "Invalid precision to cut!ep: ~D" k))  
  182.   (if (not (lisp:> (setq k (lisp:- k (bigfloat-exponent nmbr))) 0))
  183.       nmbr
  184.       (decprec! nmbr k)))
  185.  
  186. ;; This function rounds a number N at the (K+1)th place and returns an
  187. ;; equivalent number of precision K if the precision of N is greater
  188. ;; than K, else it returns the given number unchanged.
  189. (defun round!mt (nmbr k)
  190.   (unless (bigfloatp nmbr)
  191.     (error "Invalid argument to round!mt: ~S" nmbr))
  192.   (unless (and (integerp k) (not (lisp:minusp k)))
  193.     (error "Invalid precision to round!mt: ~D" k))  
  194.   (cond ((lisp:minusp (setq k (lisp:- (preci! nmbr) k 1)))
  195.      nmbr)
  196.     ((equal k 0) (round!last nmbr))
  197.     (t (round!last (decprec! nmbr k)))))
  198.  
  199. ;; This function rounds a number N and returns an equivalent number
  200. ;; having the exponent K if the exponent of N is less than K, else it
  201. ;; returns the given number unchanged.
  202. (defun round!ep (nmbr k)
  203.   (unless (bigfloatp nmbr)
  204.     (error "Invalid argument to cut!ep: ~S" nmbr))
  205.   (unless (integerp k)
  206.     (error "Invalid precision to cut!ep: ~D" k))  
  207.   (cond ((lisp:< (setq k (lisp:- (lisp:1- k) (bigfloat-exponent nmbr))) 0)
  208.      nmbr)
  209.     ((equal k 0) (round!last nmbr))
  210.     (t (round!last (decprec! nmbr k)))))
  211.  
  212. ;; This function rounds a number N at its last place.
  213. (defun round!last (nmbr)
  214.   (let ((abs-nmbr (abs (bigfloat-mantissa nmbr)))
  215.     n)
  216.     (setq n (if (lisp:< (rem abs-nmbr 10) 5)
  217.         (lisp:truncate abs-nmbr 10)
  218.         (1+ (lisp:truncate abs-nmbr 10))))
  219.     (if (lisp:minusp (bigfloat-mantissa nmbr)) (setq n (lisp:- n)))
  220.     (make-bigfloat *domain*  n (lisp:1+ (bigfloat-exponent nmbr)))))
  221.  
  222. ;;; Routines for reading/printing numbers
  223.  
  224. ;; This function reads a long number N represented by a list in a way
  225. ;; described below, and constructs a big-float representation of N.
  226. ;; Using this function, you can input any long floating-point numbers
  227. ;; without difficulty.  L is a list of integers, the first element of
  228. ;; which gives the order of N and all the next elements when
  229. ;; concatenated give the mantissa of N.
  230. ;;    **** ORDER(N)=K IF 10**K <= ABS(N) < 10**(K+1).
  231.  
  232. (defun read!lnum (l)
  233.   (loop for q in l
  234.     unless (integerp q)
  235.       do (error "Invalid argument to read!lnum: ~S" q))
  236.   (loop for term in (rest l)
  237.     for k = (lisp:ceiling (integer-length term) #.(lisp:log 10 2))
  238.     with mt = 0
  239.     and ep = (1+ (first l))
  240.     do (setq mt (lisp:+ (lisp:* mt (lisp:expt 10 k)) (abs term)))
  241.        (setq ep (lisp:- ep k))
  242.     finally (return 
  243.           (make-bigfloat *domain*
  244.             (if (lisp:plusp (second l)) mt (lisp:- mt))
  245.             ep))))
  246.  
  247. ;; This function reads a long number N represented by a list in a way
  248. ;; described below, and constructs a big-float representation of N.
  249. ;; Using this function, you can input any long floating-point numbers
  250. ;; without difficulty.  L is a list of integers, the first element of
  251. ;; which gives the order of N and all the next elements when
  252. ;; concatenated give the mantissa of N.
  253. ;;  **** ORDER(N)=K IF 10**K <= ABS(N) < 10**(K+1).       
  254.  
  255. (defun read!num (n)
  256.   (let ((exponent 0))
  257.     (multiple-value-bind (integer j)
  258.     (parse-integer n :junk-allowed t)
  259.       (unless integer
  260.     (setq integer 0))
  261.       (cond ((char= (aref n j) #\.)
  262.          (multiple-value-bind (fraction i)
  263.          (parse-integer n :start (1+ j) :junk-allowed t)
  264.            (setq integer (lisp:+ (lisp:* integer
  265.                          (lisp:expt 10 (lisp:- i j 1)))
  266.                      fraction))
  267.            (decf exponent (lisp:- i j 1)))))
  268.       (make-bigfloat *domain* integer exponent))))
  269.  
  270. ;;;  Arithmetic manipulation routines
  271.  
  272. (defun bf-abs (nmbr)
  273.   (if (lisp:> (bigfloat-mantissa nmbr) 0) nmbr
  274.       (make-bigfloat *domain*
  275.     (lisp:- (bigfloat-mantissa nmbr))
  276.     (bigfloat-exponent nmbr))))
  277.  
  278. (defun bf-minus (nmbr)
  279.   (make-bigfloat *domain*
  280.     (lisp:- (bigfloat-mantissa nmbr))
  281.     (bigfloat-exponent nmbr)))
  282.  
  283. (defun bf-plus (n1 n2)
  284.   (let ((e1 (bigfloat-exponent n1)) (e2 (bigfloat-exponent n2)))
  285.     (cond ((lisp:= e1 e2)
  286.        (make-bigfloat *domain*
  287.          (lisp:+ (bigfloat-mantissa n1) (bigfloat-mantissa n2))
  288.          e1))
  289.       ((lisp:> e1 e2)
  290.        (make-bigfloat *domain*
  291.          (lisp:+ (bigfloat-mantissa (incprec! n1 (lisp:- e1 e2)))
  292.              (bigfloat-mantissa n2))
  293.          e2))
  294.       (t (make-bigfloat *domain*
  295.            (lisp:+ (bigfloat-mantissa n1)
  296.                (bigfloat-mantissa (incprec! n2 (lisp:- e2 e1))))
  297.            e1))))) 
  298.  
  299. (defun bf-difference (n1 n2)
  300.   (let ((e1 (bigfloat-exponent n1)) (e2 (bigfloat-exponent n2)))
  301.     (cond ((lisp:= e1 e2)
  302.        (make-bigfloat *domain*
  303.          (lisp:- (bigfloat-mantissa n1) (bigfloat-mantissa n2))
  304.          e1))
  305.       ((lisp:> e1 e2)
  306.        (make-bigfloat *domain*
  307.          (lisp:- (bigfloat-mantissa (incprec! n1 (lisp:- e1 e2)))
  308.              (bigfloat-mantissa n2))
  309.          e2))
  310.       (t 
  311.        (make-bigfloat *domain*
  312.            (lisp:- (bigfloat-mantissa n1)
  313.                (bigfloat-mantissa (incprec! n2 (lisp:- e2 e1))))
  314.            e1)))))
  315.  
  316. (defun bf-times (n1 n2)
  317.   (make-bigfloat *domain*
  318.     (lisp:* (bigfloat-mantissa n1) (bigfloat-mantissa n2))
  319.     (lisp:+ (bigfloat-exponent n1) (bigfloat-exponent n2)))) 
  320.  
  321. (defun bf-quotient (n1 n2 k)
  322.   (round!mt
  323.    (with-slots (mantissa exponent) (conv!mt n1 (lisp:+ k (preci! n2) 1))
  324.      (make-bigfloat *domain*
  325.        (lisp:truncate mantissa (bigfloat-mantissa n2))
  326.        (lisp:- exponent (bigfloat-exponent n2))))
  327.    k)) 
  328.  
  329. ;; This function calculates the kth power of 'n'. The result will
  330. ;; become a long number if abs(k) >> 1.                             
  331. (defun bf-expt (number k precision)
  332.   (if (lisp:< k 0)
  333.       (/ (make-bigfloat *domain* 1 0)
  334.      (bf-expt number (lisp:- k) precision))
  335.       (%funcall (repeated-squaring
  336.               (lambda (a b) (round!mt (bf-times a b) precision))
  337.               (make-bigfloat *domain* 1 0))
  338.             number k)))
  339.  
  340. ;; This function calculates the integer quotient of 'n1' and 'n2',
  341. ;; just as the quotient" for integers does.
  342. (defun bf-floor (n1 n2)
  343.   (let ((e1 (bigfloat-exponent n1))
  344.     (e2 (bigfloat-exponent n2)))
  345.     (cond ((lisp:= e1 e2)
  346.        (make-bigfloat *domain*
  347.          (lisp:truncate (bigfloat-mantissa n1) (bigfloat-mantissa n2))
  348.          0))
  349.       ((lisp:> e1 e2)
  350.        (bf-floor (incprec! n1 (lisp:- e1 e2)) n2))
  351.       (t  (bf-floor n1 (incprec! n2 (lisp:- e2 e1)))))))
  352.  
  353. (defun bf-integer-part (num)
  354.   (with-slots (exponent mantissa) num
  355.     (if (lisp:zerop exponent) mantissa
  356.     (lisp:* mantissa (lisp:expt 10 exponent)))))
  357.  
  358. ;; This returns a lisp integer as its first return value  (perhaps
  359. ;; this should be a floating point integer)?
  360.  
  361. (defmethod floor ((number bigfloat) &optional modulus)
  362.   (let ((domain (domain-of number))
  363.     quo)
  364.     (bind-domain-context domain
  365.       (cond ((null modulus)
  366.          (setq quo (cut!ep number 0))
  367.          (values (bf-integer-part quo)
  368.              (bf-difference number quo)))
  369.         (t (setq modulus (coerce modulus (domain-of number)))
  370.            (setq quo (bf-floor number modulus))
  371.            (values (bf-integer-part quo)
  372.                (bf-difference number (bf-times quo modulus))))))))
  373.  
  374. (defmethod ceiling ((number bigfloat) &optional modulus)
  375.   (let ((domain (domain-of number))
  376.     quo)
  377.     (bind-domain-context domain
  378.       (cond ((null modulus)
  379.          (setq quo (cut!ep number 0))
  380.          (unless (eql quo number)
  381.            (setq quo (+ 1 quo)))
  382.          (values (bf-integer-part quo)
  383.              (bf-difference number quo)))
  384.         (t (setq modulus (coerce modulus (domain-of number)))
  385.            (setq quo (bf-floor number modulus))
  386.            (unless (eql quo (cut!ep quo 0))
  387.          (setq quo (+ 1 quo)))
  388.            (values (bf-integer-part quo)
  389.                (bf-difference number (bf-times quo modulus))))))))
  390.  
  391. (defmethod round ((number bigfloat) &optional modulus)
  392.   (let ((domain (domain-of number))
  393.     quo)
  394.     (bind-domain-context domain
  395.       (cond ((null modulus)
  396.          (setq quo (floor (+ (coerce 0.5 domain) number))) 
  397.          (values quo
  398.              (bf-difference number (coerce quo domain))))
  399.         (t (setq modulus (coerce modulus domain))
  400.            (setq quo (bf-floor (+ number (* (coerce 0.5 domain) modulus))
  401.                    modulus))
  402.            (values quo
  403.                (bf-difference
  404.               number
  405.               (bf-times (coerce quo domain) modulus))))))))
  406.  
  407. (defmethod truncate ((number bigfloat) &optional modulus)
  408.   (if (plus? number)
  409.       (floor number modulus)
  410.       (ceiling number modulus)))
  411.  
  412.  
  413. ;;;  Arithmetic predicates
  414.  
  415. ;; This function returns t if 'n1' > 'n2' else returns nil. 
  416. (defmethod-binary > bigfloat (n1 n2)
  417.   (with-slots ((e1 exponent)) n1
  418.     (with-slots ((e2 exponent)) n2
  419.       (cond ((lisp:=  e1 e2)
  420.          (lisp:> (bigfloat-mantissa n1) (bigfloat-mantissa n2)))
  421.         ((lisp:> e1 e2)
  422.          (lisp:> (bigfloat-mantissa (incprec! n1 (lisp:- e1 e2)))
  423.              (bigfloat-mantissa n2)))
  424.         ((lisp:> (bigfloat-mantissa n1)
  425.              (bigfloat-mantissa (incprec! n2 (lisp:- e2 e1))))
  426.          t)
  427.         (t nil)))))
  428.  
  429. (defmethod-binary < bigfloat (n1 n2)
  430.   (with-slots ((e1 exponent)) n1
  431.     (with-slots ((e2 exponent)) n2
  432.       (cond ((lisp:=  e1 e2)
  433.          (lisp:< (bigfloat-mantissa n1) (bigfloat-mantissa n2)))
  434.         ((lisp:< e1 e2)
  435.          (lisp:< (bigfloat-mantissa (incprec! n1 (lisp:- e1 e2)))
  436.              (bigfloat-mantissa n2)))
  437.         ((lisp:< (bigfloat-mantissa n1)
  438.              (bigfloat-mantissa (incprec! n2 (lisp:- e2 e1))))
  439.          t)
  440.         (t nil)))))
  441.  
  442. (defmethod-binary = bigfloat (n1 n2)
  443.   (with-slots ((e1 exponent)) n1
  444.     (with-slots ((e2 exponent)) n2
  445.       (and (lisp:=  e1 e2)
  446.        (lisp:= (bigfloat-mantissa n1) (bigfloat-mantissa n2))))))
  447.  
  448. (defmethod-binary >= bigfloat (n1 n2)
  449.   (with-slots ((e1 exponent)) n1
  450.     (with-slots ((e2 exponent)) n2
  451.       (cond ((lisp:=  e1 e2)
  452.          (lisp:>= (bigfloat-mantissa n1) (bigfloat-mantissa n2)))
  453.         ((lisp:> e1 e2)
  454.          (lisp:> (bigfloat-mantissa (incprec! n1 (lisp:- e1 e2)))
  455.              (bigfloat-mantissa n2)))
  456.         ((lisp:>= (bigfloat-mantissa n1)
  457.               (bigfloat-mantissa (incprec! n2 (lisp:- e2 e1))))
  458.          t)
  459.         (t nil)))))
  460.  
  461. (defmethod-binary <= bigfloat (n1 n2)
  462.   (with-slots ((e1 exponent)) n1
  463.     (with-slots ((e2 exponent)) n2
  464.       (cond ((lisp:=  e1 e2)
  465.          (lisp:<= (bigfloat-mantissa n1) (bigfloat-mantissa n2)))
  466.         ((lisp:< e1 e2)
  467.          (lisp:< (bigfloat-mantissa (incprec! n1 (lisp:- e1 e2)))
  468.              (bigfloat-mantissa n2)))
  469.         ((lisp:<= (bigfloat-mantissa n1)
  470.               (bigfloat-mantissa (incprec! n2 (lisp:- e2 e1))))
  471.          t)
  472.         (t nil)))))
  473.  
  474. (defun bf-integerp (x)
  475.   (and (bigfloatp x)
  476.        (not (lisp:minusp (bigfloat-exponent x)))))
  477.   
  478. (defmethod minus? ((x bigfloat))
  479.   (lisp:minusp (bigfloat-mantissa x)))
  480.  
  481. (defmethod plus? ((x bigfloat))
  482.   (lisp:plusp (bigfloat-mantissa x)))
  483.  
  484. (defmethod 0? ((number bigfloat))
  485.   (equal (bigfloat-mantissa number) 0))
  486.  
  487. (defmethod 1? ((number bigfloat))
  488.   (and (equal (bigfloat-mantissa number) 1)
  489.        (eql (bigfloat-exponent number) 0)))
  490.  
  491. (defmethod abs ((number bigfloat))
  492.   (bind-domain-context (domain-of number)
  493.     (bf-abs number)))
  494.  
  495. (defmethod minus ((number bigfloat))
  496.   (bind-domain-context (domain-of number)
  497.     (bf-minus number)))
  498.  
  499. (defmethod-binary plus bigfloat (x y)
  500.   (bind-domain-context domain
  501.     (round!mt (bf-plus x y) (fp-precision domain))))
  502.  
  503. (defmethod plus ((x bigfloat) (y float))
  504.   (plus x (coerce y (domain-of x))))
  505.  
  506. (defmethod plus ((x float) (y bigfloat))
  507.   (plus (coerce x (domain-of y)) y))
  508.  
  509. (defmethod-binary difference bigfloat (x y)
  510.   (bind-domain-context domain
  511.     (round!mt (bf-difference x y) (fp-precision domain))))
  512.  
  513. (defmethod difference ((x bigfloat) (y float))
  514.   (difference x (coerce y (domain-of x))))
  515.  
  516. (defmethod difference ((x float) (y bigfloat))
  517.   (difference (coerce x (domain-of y)) y))
  518.  
  519. (defmethod-binary times bigfloat (x y)
  520.   (bind-domain-context domain
  521.     (round!mt (bf-times x y) (fp-precision domain))))
  522.  
  523. (defmethod times ((x bigfloat) (y float))
  524.   (times x (coerce y (domain-of x))))
  525.  
  526. (defmethod times ((x float) (y bigfloat))
  527.   (times (coerce x (domain-of y)) y))
  528.  
  529. (defmethod-binary quotient bigfloat (x y)
  530.   (bind-domain-context domain
  531.     (round!mt (bf-quotient x y (fp-precision domain))
  532.           (fp-precision domain))))
  533.  
  534. (defmethod quotient ((x bigfloat) (y float))
  535.   (quotient x (coerce y (domain-of x))))
  536.  
  537. (defmethod quotient ((x float) (y bigfloat))
  538.   (quotient (coerce x (domain-of y)) y))
  539.  
  540. (defmethod expt ((number bigfloat) (k integer))
  541.   (cond ((eql k 0) (make-bigfloat (domain-of number) 1 0))
  542.     ((eql k 1) number)
  543.     (t (let ((domain (domain-of number)))
  544.          (bind-domain-context domain
  545.            (bf-expt number k (fp-precision domain)))))))
  546.  
  547.  
  548. ;; Elementary Constants
  549.  
  550. ;; This function returns the value of constant CNST of the precision
  551. ;; K, if it was calculated previously with, at least, the precision K,
  552. ;; else it returns :NOT-FOUND.
  553. (defun get!const (cnst k)
  554.   (unless (atom cnst)
  555.     (error "Invalid argument to get!const: ~S" cnst))
  556.   (unless (and (integerp k) (> k 0))
  557.     (error "Invalid precision to get!const: ~D" k))  
  558.   (let ((u (get cnst 'save!c)))
  559.     (cond ((or (null u) (< (car u) k)) nil)
  560.       ((equal (car u) k) (cdr u))
  561.       (t (round!mt (cdr u) k))))) 
  562.  
  563. ;; This function saves the value of constant CNST for the later use.
  564. (defun save!const (cnst nmbr)
  565.   (unless (atom cnst)
  566.     (error "Invalid constant for save!const: ~S" cnst))
  567.   (unless (bigfloatp nmbr)
  568.     (error "Invalid argument to  save!const: ~S" nmbr))
  569.   (setf (get cnst 'save!c) (cons (preci! nmbr) nmbr)))
  570.  
  571. ;; This function sets the value of constant CNST.  CNST is the name of
  572. ;; the constant.  L is a list of integers, which represents the value
  573. ;; of the constant in the way described in the function READ!LNUM.
  574. (defmacro set!const (constant digits)
  575.   `(progn (save!const ',constant (read!lnum ',digits))
  576.       ',constant))
  577.  
  578. (set!const !pi
  579.   (0 31415926535897932384626433832795028841971693993751058209749))
  580.  
  581. (set!const !e
  582.   (0 27182818284590452353602874713526624977572470936999595749669))
  583.  
  584. (defmacro define-bfloat-constant (name &body form)
  585.   `(defun ,name (prec)
  586.      (let ((u (get!const ',name prec)))
  587.        (when (eq u :not-found)
  588.      (setq u ,@form)
  589.      (save!const ',name u))
  590.        u)))
  591.  
  592. (defmethod pi-value ((domain real-numbers))
  593.   (bind-domain-context domain
  594.     (bf-pi (fp-precision domain))))
  595.  
  596. (defun bf-pi (precision)
  597.   (cond ((< precision 20)
  598.      (round!mt (make-bigfloat *domain* 314159265358979323846 -20)
  599.            precision))
  600.     ((get!const '!pi precision))
  601.     ((< precision 1000) (bf-pi-machin precision))
  602.     (t (bf-pi-agm precision))))
  603.  
  604. ;;  This function calculates the value of  `pi' with the precision K by       
  605. ;; using Machin's identity:           
  606. ;;          PI = 16*ATAN(1/5) - 4*ATAN(1/239).         
  607. ;; The calculation is performed mainly on integers.  
  608. (defun bf-pi-machin (k)
  609.   (let* ((k+3 (+ k 3))
  610.      s
  611.      (ss (lisp:truncate (expt 10 k+3) 5))
  612.      (n ss)
  613.      (m 1)
  614.      (x -25)
  615.      u)
  616.     (loop while (not (lisp:zerop n)) do
  617.       (setq n (lisp:truncate n x))
  618.       (setq ss (+ ss (lisp:truncate n (setq m (+ m 2))))))
  619.     (setq s (setq n (lisp:truncate (expt 10 k+3) 239)))
  620.     (setq x (- (expt 239 2)))
  621.     (setq m 1)
  622.     (loop while (not (lisp:zerop n)) do
  623.       (setq n (lisp:truncate n x))
  624.       (setq s (+ s (lisp:truncate n (setq m (+ m 2))))))
  625.     (setq u (round!mt (make-bigfloat *domain* (- (* 16 ss) (* 4 s)) (- k+3))
  626.               k))
  627.     (save!const '!pi u)
  628.     u))
  629.  
  630. ;; This function calculates the value of 'PI', with the precision K, by
  631. ;; the arithmetic-geometric mean method.  (R. Brent, JACM vol.23, #2,
  632. ;; pp.242-251(1976).)
  633. (defun bf-pi-agm (k)
  634.   (let* ((n 1)
  635.      (k2 (+ k 2))
  636.      (u (coerce 0.25 *domain*))
  637.      (half (coerce 0.5 *domain*))
  638.      (dcut (make-bigfloat *domain* 10 (- k2)))
  639.      (x (coerce 1 *domain*))
  640.      (y (bf-quotient x (bf-sqrt (coerce 2 *domain*) k2) k2))
  641.      v)
  642.     (loop while (> (bf-abs (bf-difference x y)) dcut) do
  643.       (setq v x)
  644.       (setq x (bf-times (bf-plus x y) half))
  645.       (setq y (bf-sqrt (cut!ep (bf-times y v) (- k2)) k2))
  646.       (setq v (bf-difference x v))
  647.       (setq v (bf-times (bf-times v v) (coerce n *domain*)))
  648.       (setq u (bf-difference u (cut!ep v (- k2))))
  649.       (setq n (* 2 n)))
  650.     (setq v (cut!mt (bf-expt (bf-plus x y) 2 k2) k2))
  651.     (setq u (bf-quotient v (bf-times (coerce 4 *domain*) u) k))
  652.     (save!const '!pi u)
  653.     u))
  654.  
  655. ;; This function calculates the value of 'E', the base of the natural
  656. ;; logarithm, with precision K, by summing the Taylor series for
  657. ;; EXP(X=1).
  658. (defmethod e-value ((domain real-numbers))
  659.   (bind-domain-context domain
  660.     (bf-e (fp-precision domain))))
  661.  
  662. (defun bf-e (precision)
  663.   (cond ((not (> precision 20))
  664.      (round!mt (make-bigfloat *domain* 271828182845904523536 -20)
  665.            precision))
  666.     (t (let* ((u (get!const '!e precision))
  667.           (k2 (+ precision 2))
  668.           (m 1)
  669.           (n (expt 10 k2))
  670.           (ans 0))
  671.          (when (null u)
  672.            (loop while (not (lisp:zerop n))
  673.              do (incf ans (setq n (lisp:truncate n (incf m)))))
  674.            (setq ans (+ ans (* 2 (expt 10 k2))))
  675.            (setq u (round!mt (make-bigfloat *domain* ans (- k2))
  676.                  precision))
  677.            (save!const '!e u))
  678.          u))))
  679.  
  680. ;;;  Elementary Functions. 
  681.  
  682. ;; This function calculates the square root of X with the precision K,
  683. ;; by Newton's iteration method.
  684. (defmethod sqrt ((number bigfloat))
  685.   (bind-domain-context (domain-of number)
  686.     (bf-sqrt number (fp-precision *domain*))))
  687.  
  688. (defun bf-sqrt (x k)
  689.   (if (0? x) (coerce 0 *domain*)
  690.       (let* ((k2 (+ k 2))
  691.          (ncut (- k2 (lisp:truncate (1+ (order! x)) 2)))
  692.          (half (coerce 0.5 *domain*))
  693.          (dcut (make-bigfloat *domain* 10 (- ncut)))
  694.          (dy (make-bigfloat *domain* 20 (- ncut)))
  695.          (nfig 1)
  696.          (y0 (conv!mt x 2))
  697.          y u)
  698.     (setq y0 (if (lisp:zerop (rem (bigfloat-exponent y0) 2))
  699.              (make-bigfloat *domain*
  700.                (+ 3 (* 2 (lisp:truncate (bigfloat-mantissa y0) 25)))
  701.                (lisp:truncate (bigfloat-exponent y0) 2))
  702.              (make-bigfloat *domain*
  703.                (+ 10 (* 2 (lisp:truncate (bigfloat-mantissa y0) 9)))
  704.                (lisp:truncate (- (bigfloat-exponent y0) 1) 2))))
  705.     (loop while (or (< nfig k2)
  706.             (> (bf-abs dy) dcut))
  707.           do (if (> (setq nfig (* 2 nfig)) k2)
  708.              (setq nfig k2))
  709.          (setq u (bf-quotient x y0 nfig))
  710.          (setq y (bf-times (bf-plus y0 u) half))
  711.          (setq dy (bf-difference y y0))
  712.          (setq y0 y))
  713.     (round!mt y k)))) 
  714.  
  715.  
  716. ;; This function calculates the value of the exponential function at
  717. ;; the point 'x', with the precision k, by summing terms of the Taylor
  718. ;; series for exp(z), 0 < z < 1.
  719. (defmethod exp ((number bigfloat))
  720.   (bind-domain-context (domain-of number)
  721.     (bf-exp number (fp-precision *domain*))))
  722.  
  723. (defun bf-exp (x k)
  724.   (cond ((0? x) (coerce 1 *domain*))
  725.     (t (let* ((k2 (+ k 2))
  726.           (one (coerce 1 *domain*))
  727.           (y (bf-abs x))
  728.           (m (floor y))
  729.           (q (coerce m *domain*))
  730.           (r (bf-difference y q))
  731.           yq yr)
  732.          (setq yq (if (lisp:zerop m) one
  733.               (bf-expt (bf-e k2) m k2)))
  734.          (cond ((0? r) (setq yr one))
  735.            (t (let ((j 0) (n 0)
  736.                 (dcut (make-bigfloat *domain* 10 (- k2)))
  737.                 (ri one) (tm one)
  738.                 fctrial)
  739.             (setq yr one)
  740.             (setq m 1)
  741.             (loop while (> tm dcut) do
  742.               (setq fctrial
  743.                 (coerce
  744.                   (setq m (* m (setq j (1+ j)))) *domain*))
  745.               (setq ri (cut!ep (bf-times ri r) (- k2)))
  746.               (setq n (max 1 (+ (- k2 (order! fctrial))
  747.                         (order! ri))))
  748.               (setq tm (bf-quotient ri fctrial n))
  749.               (setq yr (bf-plus yr tm))
  750.               (cond ((lisp:zerop (rem j 10))
  751.                  (setq yr (cut!ep yr (- k2)))))))))
  752.          (setq y (cut!mt (bf-times yq yr) (1+ k)))
  753.          (if (minus? x) (bf-quotient one y k)
  754.          (round!last y)))))) 
  755.  
  756.  
  757. ;; This function calculates log(x) by summing terms of the     
  758. ;; Taylor series for LOG(1+Z), 0 < Z < 0.10518. 
  759. (defmethod log ((x bigfloat))
  760.   (let ((domain (domain-of x)))
  761.     (bind-domain-context domain
  762.       (bf-log x (fp-precision domain)))))
  763.  
  764. (defmethod-binary log2 bigfloat (x base)
  765.   (let ((k2 (+ 2 (fp-precision domain))))
  766.     (bind-domain-context domain
  767.       (bf-quotient (bf-log x k2) (bf-log base k2) (- k2 2)))))
  768.  
  769. (defun bf-log (x k)
  770.   (when (not (plus? x))
  771.     (error "Invalid argument to log: ~S" x))
  772.   (unless (and (integerp k) (> k 0))
  773.     (error "Invalid precision to log: ~D" k))
  774.   (cond ((= x 1)
  775.      (coerce 0 *domain*))
  776.     (t (let* ((m 0)
  777.           (k2 (+ k 2))
  778.           (one (coerce 1 *domain*))
  779.           (ee (bf-e k2))
  780.           (es (bf-exp (coerce 0.1 *domain*) k2))
  781.           sign l y z)
  782.          (cond ((> x one) (setq sign one) (setq y x))
  783.            (t (setq sign (bf-minus one))
  784.               (setq y (bf-quotient one x k2))))
  785.          (cond ((< y ee)
  786.             (setq z y))
  787.            (t 
  788.             (cond ((lisp:zerop (setq m (lisp:truncate (* (order! y) 23)
  789.                                   10)))
  790.                (setq z y))
  791.               (t (setq z (bf-quotient y (bf-expt ee m k2) k2))))
  792.             (loop while (> z ee) do
  793.               (setq m (1+ m))
  794.               (setq z (bf-quotient z ee k2)))))
  795.          (setq l (coerce m *domain*))
  796.          (setq y (coerce 0.1 *domain*))
  797.          (loop while (> z es) do
  798.            (setq l (bf-plus l y))
  799.            (setq z (bf-quotient z es k2)))
  800.          (setq z (bf-difference z one))
  801.          (prog (n dcut tm zi)
  802.            (setq n 0)
  803.            (setq y (setq tm (setq zi z)))
  804.            (setq z (bf-minus z))
  805.            (setq dcut (make-bigfloat *domain* 10 (- k2)))
  806.            (setq m 1)
  807.            (loop while (> (bf-abs tm) dcut) do
  808.          (setq zi (cut!ep (bf-times zi z) (- k2)))
  809.          (setq n (max 1 (+ k2 (order! zi))))
  810.          (setq tm
  811.                (bf-quotient zi (coerce (setq m (1+ m)) *domain*)
  812.               n))
  813.          (setq y (bf-plus y tm))
  814.          (cond ((lisp:zerop (rem m 10))
  815.             (setq y (cut!ep y (- k2)))))))
  816.          (setq y (bf-plus y l))
  817.          (round!mt (bf-times sign y) k)))))
  818.  
  819. ;; This function calculates log(x), the value of the logarithmic
  820. ;; function at the point 'x', with the precision k, by solving x =
  821. ;; exp(y) by Newton's method.
  822. ;;  x > 0, k is a positive integer                        
  823. #+Ignore
  824. (defun bf-log-newton (x k)
  825.   (when (not (plus? x))
  826.     (error "Invalid argument to log: ~S" x))
  827.   (unless (and (integerp k) (> k 0))
  828.     (error "Invalid precision to log: ~D" k))
  829.   (cond ((bf-equal x (coerce 1 *domain*))
  830.      (coerce 0 *domain*))
  831.     (t (let* ((k2 (+ k 2))
  832.           m (one (coerce 1 *domain*))
  833.           (ee (bf-e (+ k2 2)))
  834.           sign y z)
  835.          (cond ((> x one)
  836.             (setq sign one)
  837.             (setq y x))
  838.            (t (setq sign (bf-minus one))
  839.               (setq y (bf-quotient one x k2))))
  840.          (if (< y ee)
  841.          (setq m 0 z y)
  842.          (cond ((lisp:zerop
  843.              (setq m (lisp:truncate (lisp:* (order! y) 23) 10)))
  844.             (setq z y))
  845.                (t (setq z (bf-quotient y (bf-expt ee m k2) k2))
  846.               (loop while (> z ee) do
  847.                 (setq m (1+ m))
  848.                 (setq z (bf-quotient z ee k2))))))
  849.          (let ((nfig 0) (n 0)
  850.            (dcut (make-bigfloat *domain* 10 (- k2)))
  851.            dx
  852.            (dy (make-bigfloat *domain* 20 (- k2)))
  853.            x0)
  854.            (setq y (bf-quotient (bf-difference z one)
  855.                     (coerce 1.72 *domain*) 2))
  856.            (setq nfig 1)
  857.            (loop while (or (< nfig k2) (> (bf-abs dy) dcut)) do
  858.          (cond
  859.            ((> (setq nfig (* 2 nfig)) k2)
  860.             (setq nfig k2)))
  861.          (setq x0 (exp* y nfig))
  862.          (setq dx (bf-difference z x0))
  863.          (setq n (max 1 (+ nfig (order! dx))))
  864.          (setq dy (bf-quotient dx x0 n))
  865.          (setq y (bf-plus y dy))))
  866.          (setq y (bf-plus (coerce m *domain*) y))
  867.          (round!mt (bf-times sign y) k)))))
  868.  
  869. ;; This function calculates sin(x), the value of the sine function at
  870. ;; the point 'x', with the precision k, by summing terms of the Taylor
  871. ;; series for:   sin(z), 0 < Z < PI/4.    
  872. (defmethod sin ((number bigfloat))
  873.   (let ((domain (domain-of number)))
  874.     (bind-domain-context domain
  875.       (bf-sin number (fp-precision domain)))))
  876.  
  877. (defun bf-sin (x k)
  878.   (cond ((0? x) (coerce 0 *domain*))
  879.     ((minus? x) (bf-minus (bf-sin (bf-minus x) k)))
  880.     (t (let* ((k2 (+ k 2))
  881.           (m (preci! x))
  882.           (pi4 (bf-times (bf-pi (+ k2 m)) (coerce 0.25 *domain*)))
  883.           sign q r y)
  884.          (cond ((< x pi4)
  885.             (setq m 0)
  886.             (setq r x))
  887.            (t (setq m (floor (setq q (bf-floor x pi4))))
  888.               (setq r (bf-difference x (bf-times q pi4)))))
  889.          (setq sign (coerce 1 *domain*))
  890.          (cond ((not (< m 8)) (setq m (rem m 8))))
  891.          (cond ((not (< m 4))
  892.             (setq sign (bf-minus sign))
  893.             (setq m (- m 4))))
  894.          (cond ((equal m 1)
  895.             (setq r (cut!mt (bf-difference pi4 r) k2))
  896.             (bf-times sign (bf-cos r k)))
  897.            ((equal m 2)
  898.             (setq r (cut!mt r k2))
  899.             (bf-times sign (bf-cos r k)))
  900.            (t (unless (equal m 0)
  901.             (setq r (cut!mt (bf-difference pi4 r) k2)))
  902.               (let* ((ncut (- k2 (min 0 (1+ (order! r)))))
  903.                  (dcut (make-bigfloat *domain* 10 (- ncut)))
  904.                  (tm r) (ri r) (j 1) n fctrial)
  905.             (setq y r)
  906.             (setq r (bf-minus (cut!ep (bf-times r r) (- ncut))))
  907.             (setq m 1)
  908.             (loop while (> (bf-abs tm) dcut) do
  909.               (setq j (+ j 2))
  910.               (setq fctrial
  911.                 (coerce
  912.                   (setq m (* m j (1- j))) *domain*))
  913.               (setq ri (cut!ep (bf-times ri r) (- ncut)))
  914.               (setq n (max 1 (+ (- k2 (order! fctrial)) (order! ri))))
  915.               (setq tm (bf-quotient ri fctrial n))
  916.               (setq y (bf-plus y tm))
  917.               (cond ((lisp:zerop (rem j 20))
  918.                  (setq y (cut!ep y (- ncut)))))))
  919.               (round!mt (bf-times sign y) k)))))))
  920.  
  921. ;; This function calculates cos(x), the value of the cosine function at
  922. ;; the point 'x', with the precision k, by summing terms of the Taylor
  923. ;; series for:   cos(z), 0 < Z < PI/4.    
  924. (defmethod cos ((number bigfloat))
  925.   (let ((domain (domain-of number)))
  926.     (bind-domain-context domain
  927.       (bf-cos number (fp-precision domain)))))
  928.  
  929. (defun bf-cos (x k)
  930.   (unless (and (integerp k) (> k 0))
  931.     (error "Invalid precision to cos: ~D" k))
  932.   (cond ((0? x) (coerce 1 *domain*))
  933.     (t (when (minus? x)
  934.          (setq x  (- x)))
  935.        (let* ((k2 (+ k 2))
  936.           (m (preci! x))
  937.           (pi4 (/ (bf-pi (+ k2 m)) 4))
  938.           sign q r y)
  939.          (cond ((< x pi4)
  940.             (setq m 0)
  941.             (setq r x))
  942.            (t (setq m (floor (setq q (floor x pi4))))
  943.               (setq r (- x (* q pi4)))))
  944.          (setq sign (coerce 1 *domain*))
  945.          (cond ((not (< m 8)) (setq m (rem m 8))))
  946.          (cond ((not (< m 4))
  947.             (setq sign (- sign))
  948.             (setq m (- m 4))))
  949.          (cond ((not (< m 2)) (setq sign (- sign))))
  950.          (cond ((equal m 1)
  951.             (setq r (cut!mt (- pi4 r) k2))
  952.             (bf-times sign (bf-sin r k)))
  953.            ((equal m 2)
  954.             (setq r (cut!mt r k2))
  955.             (bf-times sign (bf-sin r k)))
  956.            (t (when (= m 3)
  957.             (setq r (cut!mt (- pi4 r) k2)))
  958.               (let ((j 0) (n 0)
  959.                 (dcut (make-bigfloat *domain* 10 (- k2)))
  960.                 fctrial ri tm)
  961.             (setq y (setq ri (setq tm (coerce 1 *domain*))))
  962.             (setq r (- (cut!ep (* r r) (- k2))))
  963.             (setq m 1)
  964.             (loop while (> (bf-abs tm) dcut) do
  965.               (setq j (+ j 2))
  966.               (setq fctrial
  967.                 (coerce
  968.                   (setq m (* m j (- j 1))) *domain*))
  969.               (setq ri (cut!ep (* ri r) (- k2)))
  970.               (setq n (max 1 (+ (- k2 (order! fctrial))
  971.                         (order! ri))))
  972.               (setq tm (bf-quotient ri fctrial n))
  973.               (setq y (+ y tm))
  974.               (cond ((equal (rem j 20) 0)
  975.                  (setq y (cut!ep y (- k2)))))))
  976.               (round!mt (* sign y) k)))))))
  977.  
  978. ;; This function calculates tan(x), the value of the tangent function
  979. ;; at the point 'x', with the precision k, by calculating       
  980. ;;          sin(x)  or  cos(x) = sin(pi/2-x).       
  981. (defmethod tan ((number bigfloat))
  982.   (let ((domain (domain-of number)))
  983.     (bind-domain-context domain
  984.       (bf-tan number (fp-precision domain)))))
  985.  
  986. (defun bf-tan (x k)
  987.   (unless (and (integerp k) (> k 0))
  988.     (error "Invalid precision to tan: ~D" k))
  989.   (cond ((0? x) (coerce 0 *domain*))
  990.     ((minus? x) (bf-minus (bf-tan (bf-minus x) k)))
  991.     (t (let* ((k2 (+ k 2))
  992.           (one (coerce 1 *domain*))
  993.           (m (preci! x))
  994.           (pi4 (bf-times (bf-pi (+ k2 m)) (coerce 0.25 *domain*)))
  995.           sign q r)
  996.          (cond ((< x pi4)
  997.             (setq m 0)
  998.             (setq r x))
  999.            (t (setq m (floor (setq q (bf-floor x pi4))))
  1000.               (setq r (bf-difference x (bf-times q pi4)))))
  1001.          (cond ((not (< m 4)) (setq m (rem m 4))))
  1002.          (setq sign (if (< m 2) one (bf-minus one)))
  1003.          (cond ((or (= m 1) (= m 3))
  1004.             (setq r (bf-difference pi4 r))))
  1005.          (setq r (cut!mt r k2))
  1006.          (cond ((or (equal m 0) (equal m 3))
  1007.             (setq r (bf-sin r k2))
  1008.             (setq q (bf-difference one (bf-times r r)))
  1009.             (setq q (bf-sqrt (cut!mt q k2) k2))
  1010.             (bf-times sign (bf-quotient r q k)))
  1011.            (t (setq r (bf-sin r k2))
  1012.               (setq q (bf-difference one (bf-times r r)))
  1013.               (setq q (bf-sqrt (cut!mt q k2) k2))
  1014.               (bf-times sign (bf-quotient q r k))))))))
  1015.  
  1016. ;; This function calculates asin(x), the value of the arcsine function
  1017. ;; at the point 'x', with the precision k, by calculating        
  1018. ;;          atan(x/sqrt(1-x**2))  
  1019. ;; The answer is in the range <-pi/2 , pi/2>.  
  1020. (defmethod asin ((number bigfloat))
  1021.   (let ((domain (domain-of number)))
  1022.     (bind-domain-context domain
  1023.       (bf-asin number (fp-precision domain)))))
  1024.  
  1025. (defun bf-asin (x k)
  1026.   (when (or (> (bf-abs x) (coerce 1 *domain*)))
  1027.     (error "Invalid argument to asin: ~S" x))
  1028.   (unless (and (integerp k) (> k 0))
  1029.     (error "Invalid precision to asin: ~D" k))
  1030.   (cond ((minus? x) (bf-minus (bf-asin (bf-minus x) k)))
  1031.     (t (let ((k2 (+ k 2))
  1032.          (one (coerce 1 *domain*)))
  1033.          (cond ((< (bf-difference one x)
  1034.                (make-bigfloat *domain* 10 (- k2)))
  1035.             (round!mt (bf-times (bf-pi (1+ k)) (coerce 0.5 *domain*))
  1036.                   k))
  1037.            (t (bf-atan
  1038.                (bf-quotient x (bf-sqrt (cut!mt (- 1 (* x x)) k2) k2)
  1039.                     k2)
  1040.                k))))))) 
  1041.  
  1042.  
  1043. ;; This function calculates acos(x), the value of the arccosine
  1044. ;; function at the point 'x', with the precision k, by calculating        
  1045. ;;          atan(sqrt(1-x**2)/x)  if  x > 0  or      
  1046. ;;          atan(sqrt(1-x**2)/x) + pi  if  x < 0.    
  1047. ;; the answer is in the range [0 , pi).        
  1048. (defmethod acos ((number bigfloat))
  1049.   (let ((domain (domain-of number)))
  1050.     (bind-domain-context domain
  1051.       (bf-acos number (fp-precision domain)))))
  1052.  
  1053. (defun bf-acos (x k)
  1054.   (when (or (> (bf-abs x) (coerce 1 *domain*)))
  1055.     (error "Invalid argument to acos: ~S" x))
  1056.   (unless (and (integerp k) (> k 0))
  1057.     (error "Invalid precision to acos: ~D" k))
  1058.   (let ((k2 (+ k 2))
  1059.     y)
  1060.     (cond ((< (bf-abs x) (make-bigfloat *domain* 50 (- k2)))
  1061.        (round!mt (bf-times (bf-pi (+ k 1)) (coerce 0.5 *domain*))
  1062.              k))
  1063.       (t (setq y (bf-quotient
  1064.               (bf-sqrt (cut!mt
  1065.                 (bf-difference (coerce 1 *domain*)
  1066.                            (bf-times x x))
  1067.                 k2) k2)
  1068.               (bf-abs x)
  1069.               k2))
  1070.          (if (minus? x)
  1071.          (round!mt (bf-difference (bf-pi (+ k 1)) (bf-atan y k))
  1072.                k)
  1073.          (bf-atan y k))))))
  1074.  
  1075. ;; this function calculates atan(x), the value of the arctangent
  1076. ;; function at the point 'x', with the precision k, by summing terms
  1077. ;; of the Taylor series for atan(z)  if  0 < z < 0.42.  
  1078. ;;   otherwise the following identities are used!  
  1079. ;;       atan(x) = pi/2 - atan(1/x)  if  1 < x  and 
  1080. ;;       atan(x) = 2*atan(x/(1+sqrt(1+x**2)))       
  1081. ;;             if  0.42 <= x <= 1.                     
  1082. ;; the answer is in the range [-pi/2, pi/2).    
  1083. (defmethod atan ((number bigfloat) &optional base)
  1084.   (when base
  1085.     (error "Two argument atan not implemented yet"))
  1086.   (let ((domain (domain-of number)))
  1087.     (bind-domain-context domain
  1088.       (bf-atan number (fp-precision domain)))))
  1089.  
  1090. (defun bf-atan (x k)
  1091.   (unless (and (integerp k) (> k 0))
  1092.     (error "Invalid precision to atan: ~D" k))
  1093.   (cond ((0? x) (coerce 0 *domain*))
  1094.     ((minus? x) (bf-minus (bf-atan (bf-minus x) k)))
  1095.     (t (let* ((k2 (+ k 2))
  1096.           (one (coerce 1 *domain*))
  1097.           (pi4 (bf-times (bf-pi k2) (coerce 0.25 *domain*)))
  1098.           y z)
  1099.          (cond ((= x 1)
  1100.             (round!mt pi4 k))
  1101.            ((> x one)
  1102.             (round!mt
  1103.              (bf-difference (bf-plus pi4 pi4)
  1104.                     (bf-atan (bf-quotient one x k2) (+ k 1)))
  1105.              k))
  1106.            ((< x (coerce 0.42 *domain*))
  1107.             (let* ((m 1) (n 0)
  1108.                (ncut (- k2 (min 0 (+ (order! x) 1))))
  1109.                (dcut (make-bigfloat *domain* 10 (- ncut)))
  1110.                (zi x)
  1111.                (tm x))
  1112.               (setq y tm)
  1113.               (setq z (bf-minus (cut!ep (bf-times x x) (- ncut))))
  1114.               (loop while (> (bf-abs tm) dcut) do
  1115.             (setq zi (cut!ep (bf-times zi z) (- ncut)))
  1116.             (setq n (max 1 (+ k2 (order! zi))))
  1117.             (setq tm (bf-quotient
  1118.                     zi (coerce (setq m (+ m 2)) *domain*)
  1119.                     n))
  1120.             (setq y (bf-plus y tm))
  1121.             (cond ((lisp:zerop (rem m 20))
  1122.                    (setq y (cut!ep y (- ncut)))))))
  1123.             (round!mt y k))
  1124.            (t (setq y (bf-plus one (cut!mt (bf-times x x) k2)))
  1125.               (setq y (bf-plus one (bf-sqrt y k2)))
  1126.               (setq y (bf-atan (bf-quotient x y k2) (+ k 1)))
  1127.               (round!mt (bf-times y (coerce 2 *domain*)) k)))))))
  1128.  
  1129. ;; this function calculates arcsin(x), the value of the arcsine
  1130. ;; function at the point 'x', with the precision k, by solving
  1131. ;;    x = sin(y)  if  0 < x <= 0.72,  or       
  1132. ;;    sqrt(1-x**2) = sin(y)  if  0.72 < x,     
  1133. ;; by Newton's iteration method.               
  1134. ;; the answer is in the range [-pi/2, pi/2).
  1135. #+Ignore
  1136. (defun bf-asin-newton (x k)
  1137.   (when (or (> (bf-abs x) (coerce 1 *domain*)))
  1138.     (error "Invalid argument to asin: ~S" x))
  1139.   (unless (and (integerp k) (> k 0))
  1140.     (error "Invalid precision to asin: ~D" k))
  1141.   (cond ((0? x) (coerce 0 *domain*))
  1142.     ((minus? x) (bf-minus (bf-asin-newton (bf-minus x) k)))
  1143.     (t (let* ((k2 (+ k 2))
  1144.           (dcut (make-bigfloat *domain* 10 (+ (- k2) (order! x) 1)))
  1145.           (one (coerce 1 *domain*))
  1146.           (pi2 (bf-times (bf-pi (+ k2 2)) (coerce 0.5 *domain*)))
  1147.           y)           
  1148.          (cond ((< (- 1 x) dcut)
  1149.             (round!mt pi2 k))
  1150.            ((> x (coerce 0.72 *domain*))
  1151.             (setq y (cut!mt (bf-difference one (bf-times x x)) k2))
  1152.             (setq y (bf-asin-newton (bf-sqrt y k2) k))
  1153.             (round!mt (bf-difference pi2 y) k))
  1154.            (t (let ((nfig 1)
  1155.                 (n 0)
  1156.                 (dy one)
  1157.                 cx dx x0)
  1158.             (setq y x)
  1159.             (loop while (or (< nfig k2)
  1160.                     (> (bf-abs dy) dcut))
  1161.                   do (cond ((> (setq nfig (* 2 nfig)) k2)
  1162.                     (setq nfig k2)))
  1163.                  (setq x0 (bf-sin y nfig))
  1164.                  (setq cx (bf-sqrt (cut!mt (- 1 (* x0 x0))
  1165.                                nfig)
  1166.                          nfig))
  1167.                  (setq dx (- x x0))
  1168.                  (setq n (max 1 (+ nfig (order! dx))))
  1169.                  (setq dy (bf-quotient dx cx n))
  1170.                  (setq y (bf-plus y dy)))
  1171.             (round!mt y k))))))))
  1172.  
  1173. ;; This function calculates arccos(x), the value of the arccosine
  1174. ;; function at the point 'x', with the precision k, by calculating
  1175. ;;    arcsin(sqrt(1-x**2))  if  x > 0.72  and    
  1176. ;;    pi/2 - arcsin(x)  otherwise.
  1177. ;; The answer is in the range [0, pi).          
  1178. #+ignore
  1179. (defun bf-acos-newton (x k)
  1180.   (when (or (> (bf-abs x) (coerce 1 *domain*)))
  1181.     (error "Invalid argument to acos: ~S" x))
  1182.   (unless (and (integerp k) (> k 0))
  1183.     (error "Invalid precision to acos: ~D" k))
  1184.   (cond ((bf-<= x (coerce 0.72 *domain*))
  1185.          (round!mt
  1186.       (bf-difference
  1187.        (bf-times (bf-pi (+ k 1)) (coerce 0.5 *domain*))
  1188.        (bf-asin-newton x k))
  1189.       k))
  1190.     (t (bf-asin-newton
  1191.         (bf-sqrt
  1192.          (cut!mt (bf-difference (coerce 1 *domain*) (* x x))
  1193.              (+ k 2))
  1194.          (+ k 2))
  1195.         k))))
  1196.  
  1197. ;; This function calculates arctan(x), the value of the arctangent
  1198. ;; function at the point 'x', with the precision k, by calculating
  1199. ;;     arcsin(x/sqrt(1+x**2))
  1200. ;; The answer is in the range [-pi/2, pi/2).  
  1201. #+Ignore
  1202. (defun bf-atan-newton (x k)
  1203.   (unless (and (integerp k) (> k 0))
  1204.     (error "Invalid precision to atan: ~D" k))
  1205.   (cond ((minus? x) (bf-minus (bf-atan-newton (bf-minus x) k)))
  1206.     (t (bf-asin-newton
  1207.         (bf-quotient x (bf-sqrt (cut!mt (+ 1 (* x x)) (+ k 2))
  1208.                     (+ k 2))
  1209.              (+ k 2))
  1210.         k))))
  1211.  
  1212. (defmethod-binary expt bigfloat (x y)
  1213.   (bind-domain-context domain
  1214.     (cond ((bf-integerp y) (expt x (floor y)))
  1215.       ((minus? y)
  1216.        (/ 1 (expt x (bf-minus y))))
  1217.       (t (let ((n (fp-precision domain))
  1218.            (xp (abs x))
  1219.            yp) 
  1220.            (cond ((bf-integerp (bf-times y (coerce 2 domain)))
  1221.               (setq xp (incprec! xp 1)) 
  1222.               (setq yp (round!mt
  1223.                 (bf-times (expt xp (floor y))
  1224.                       (bf-sqrt xp (+ n 1)))
  1225.                 n)))
  1226.              (t (setq yp (bf-exp (* y (bf-log xp (1+ n))) n))))
  1227.            (cond ((minus? x) (bf-minus yp)) (t yp)))))))
  1228.  
  1229.